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ABSTRACT 



We have apphed axisymmetric magnetohydrodynamic (MHD) simulations in 
order to investigate the impact of the accretion disk magnetic flux profile on 
the collimation of jets. Using the ZEUS-3D code modified for magnetic diffusiv- 
ity, our simulations evolve from an initial state in hydrostatic equilibrium and 
\^ ', a force-free magnetic field configuration. Considering a power law for the disk 

poloidal magnetic field profile Bp ~ r~'^ and for the density profile of the disk 
l/^ i wind p ~ r^f^f, we have performed a systematic parameter study over a wide 

range of parameters fi and fip. We apply a toy parameterization for the magnetic 
diffusivity derived from the internal turbulent Alfvenic pressure. We find that 
the degree of collimation (quantified by the ratio of mass flow rates in axial and 
lateral direction) decreases for a steeper disk magnetic field profile (increasing 
^1 /i). Varying the total magnetic flux does not change the degree of jet collimation 

substantially, it only affects the time scale of outflow evolution and the terminal 
I jet speed. As our major result we find a general relation between the collimation 

^ I degree with the disk wind magnetization power law exponent. Outflows with high 

degree of collimation resulting from a flat disk magnetic field profile tend to be 
unsteady, producing axially propagating knots as discussed earlier in the litera- 
ture. Depending slightly on the inflow density profile this unsteady behavior sets 
in for p < 0.4. We also performed simulations of jet formation with artificially 
enhanced decay of the toroidal magnetic field component in order to investigate 
the idea of a purely "poloidal collimation" previously discussed in the literature. 
These outflows remain only weakly collimated and propagate with lower velocity. 
Thanks to our large numerical grid size (about 7 x 14 AU for protostars), we may 
apply our results to recently observed hints of jet rotation (DC Tau) indicating 



^Part of the numerical work has been accomphshed at the Astrophysikahsches Institut Potsdam and at 
the University of Potsdam, Germany 
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a relatively flat disk magnetic field profile, ~ 0.5. In general, our results are 
applicable to both stellar and extragalactic sources of MHD jets. 

Subject headings: accretion, accretion disks — MHD — ISM: jets and outflows 
— stars: mass loss — stars: pre-main sequence — galaxies: jets 

1. Introduction 

Astrophysical jets as highly coUimated high speed beams of matter have been observed 
as common phenomenon in a variety of astronomical sources, among them young stars, 
micro-quasars and active galactic nuclei. The current understanding of jet formation is that 
jets arc launched by magnetohydrodynamic (MHD) processes in the vicinity of the central 
jet object - an accretion disk surrounding a protostar or a black hole (Blandford & Payne 
1982; Pudritz & Norman 1983; Camenzind 1990). The general properties of jet sources and 
the ongoing physical processes have been investigated for a long time, however, the principal 
mechanism which actually launches a jet from the disk at a certain time is not yet known. 

A number of numerical simulations investigating the MHD jet formation have been pub- 
hshed. We may distinguish between those simulations which take into account the evolution 
of the disk structure and those which consider the disk as a fixed boundary condition for the 
simulation of the disk wind. The first approach allows to study directly the launch^ of the 
outfiow, i.e the ejection of accreting plasma into an outfiow in direction vertical to the disk. 
However, due to the limited time resolution and spatial coverage typical for these simula- 
tions, little can be learned concerning the ultimate jet acceleration and collimation which 
takes place at radii beyond the Alfven surface (Uchida & Shibata 1985; Miller & Stone 1997; 
Kudoh ct al. 1998). The second approach allows to investigate the long term, evolution of 
jet formation on a large spatial scale (Ouyed & Pudritz 1997a,b; Krasnopolsky et al. 1999; 
Fendt & Elstner 2000; Fendt & Cemeljic 2002). The mechanism which launches the outfiow 
out of the accretion stream cannot be studied by this approach. 

In this paper we study how jet formation - i.e. jet collimation and acceleration - depends 
on the magnetic field profile of the jet launching accretion disk. The initial magnetic field is 
located in an initially hydrostatic density distribution. We prescribe various magnetic field 
profiles along the disk surface which is taken as a boundary condition fixed in time. The 
coronal (i.e. jet) magnetic field evolves in time from the initial state and governed by a 



^By jet formation and jet launching we denote the acceleration and collimation of a disk wind, and the 
generation of a disk wind out of the accretion flow, respectively 
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prescribed mass inflow from the disk surface into the corona. We apply the ZEUS-3D MHD 
code modified for magnetic diffusivity (Fendt & Cemeljic 2002). It is clear that the "real" 
disk-jet magnetic field and the mass flux into the jet is a result of dissipative MHD processes 
in the disk in interaction with the external field structure. The disk structure and evolution 
will also be affected by the existence of a magnetized disk wind/jet. Although most of the 
processes involved are understood to some extent, the self-consistent dynamical modeling 
of such systems has just started (e.g. simulations of a disk dynamo generated jet magnetic 
magnetic field (von Rekowski & Brandenburg 2004)) and the results have to be considered 
still as preliminary. For the purpose of the present paper we have to keep in mind that some 
of the parameter runs considered may be more realistic than others concerning the internal 
disk physics. 

We note that in a contemporaneous and independent study Pudritz et al. (2006a) have 
undertaken a similar approach. Our results generally agree with their paper, but also tran- 
scend their approach in several aspects. Among the additional features treated in our paper 
are a much larger grid, the description of a turbulent magnetic diffusivity, a quantitative 
measure of the coUimation degree, and a broader range of parameters studied. 

Our paper is organized as follows. Section 2 briefiy summarizes the issue of MHD jet 
coUimation. Section 3 introduces our model approach. In Sect. 4 we present our results before 
concluding the paper with the summary. Preliminary results of the present investigation have 
been published earlier (see Fig. 2 in Riidiger (2002)). 

2. Self-collimation of rotating MHD outflows 

Theoretical studies considering analytical solutions of the stationary, axisymmetric 
MHD force-balance in the asymptotic jet flow (i.e. far from the source) have clearly shown 
that such outflows must collimate into a narrow beam when carrying a net poloidal electric 
current (Heyvaerts & Norman 1989; Chiueh et al. 1991). This fundamental statement has 
been recently conflrmed and generalized (Okamoto 2003; Heyvaerts & Norman 2003). 

2.1. Stationeiry studies of MHD jets 

A self-consistent treatment of the coUimation process - i.e. the question how coUimation 
is achieved along the initial outfiow - has to consider at least a 2.5-dimensional (i.e. ax- 
isymmetric three-dimensional) MHD problem. Until about a decade ago, MHD jet theory 
has been mostly limited to the stationary approach. Time-independent studies, however. 
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cannot really prove whether a MHD solution derived will actually be realized during the 
time evolution of a forming jet. Hence, the final answer concerning jet self-collimation can 
only be given by time- dependent MHD simulations. Although stationary self- similar models 
gave clear indication for collimated MHD jets (Blandford & Payne 1982; Sauty & Tsinganos 
1994; Li 1993; Contopoulos & Lovelace 1994) one has to keep in mind that such a constraint 
has implications which are critical for coUimation, but are not really feasible for jets. For 
example, self-similarity considers an infinite jet radius and does exclude the symmetry axis 
of the outflow. Truly 2.5-dimensional (axisymmetric 3D) solutions considering the local 
force-balance and global boundary conditions have been obtained only in the force-free limit 
(Fendt et al. 1995; Lery et al. 1998; Fendt & Memola 2001), and provide the shape of the col- 
limating jet as determined by the internal force-equilibrium. So far, the only self-consistent 
two-dimensional stationary MHD solution of a collimating wind flow has been published 
decades ago (Sakurai 1985, 1987). These solutions coUimate on logarithmic scales only - a 
result of the low rotation rate applied. 



2.2. Time-dependent MHD simulations 

The first numerical evidence of the MHD jet self-collimation process has been provided 
by Ouyed & Pudritz (1997a) following a pioneering approach for disk winds by Ustyugova et 
al. (1995). The simulations by Ouyed & Pudritz were been performed from an initially force- 
free setup and fixed boundary conditions for the mass inflow into the jet from an underlying 
Keplerian disk surface. After a few hundred of inner disk rotations the MHD disk wind 
coUimates into a narrow beam and saturates into a stationary state outflow. For a recent 
review on numerical progress in simulating disk winds and jets we refer to (Pudritz et al. 
2006b). 

It is clear that different inflow conditions for the disk wind (mass flow rate or magnetic 
field profile) will result in a different dynamical evolution of the jet, as there are different time 
scales, velocities or a different degree of coUimation. It is therefore interesting to investigate 
a wide parameter range for the leading jet parameters. This has partly been done in a 

number of papers discussing the variation of the jet mass load (Ouyed & Pudritz 1999; Kato 
et al. 2002), a time-dependent disk magnetic field inclination (Krasnopolsky et al. 1999), 
magnetic diffusivity in the jet (Kuwabara et al. 2000; Fendt & Cemeljic 2002; Kuwabara et 
al. 2005), time-dependent mass loading (Vitorino et al. 2003), the extension of the simulation 
box on coUimation (Ustyugova et al. 1999), or, in difference to these cases of a monotonously 
distributed disk magnetic field, the long-term evolution of a stellar dipolar magnetosphere 
in connection with a Keplerian disk (Fendt & Elstner 2000). 
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The next step will be to include the disk structure in the numerical simulations (Hayashi 
et al. 1996; Miller & Stone 1997; Cassc & Keppens 2002; Kudoh et al. 2002) and follow the 
jet launching process over many hundreds of disk rotations. Substantial progress have been 
made in this respect in simulations already in particular concerning the accretion process 
from disk down to the star (Romanova et al. 2002, 2003, 2004, 2005) but coUimated disk jets 
have not yet been found. 

2.3. "Toroidal or poloidal collimation"? 

Concerning jet sclf-collimation it is interesting to mention a proposal by Spruit et al. 
(1997) pointing out that the toroidally dominated jet magnetic field is kink-instahle. Spruit 
et al. suggested that jets should be rather collimatcd by the poloidal disk magnetic field 
pressure and less by the toroidal pinching force. The authors find that ''poloidal collima- 
tion" works best for a disk magnetic field profile |i?p| ~ with jj, < 1.3. Indication for 
"poloidal collimation" of outflows has been reported also by MHD simulations (Matt et al. 
2003), launching the outflow in a dipolar stellar magnetosphere and coUimating it by the 
surrounding disk magnetic fleld. 

The kink instability of jets can be investigated by non-axisymmetric simulations only. It 
is therefore important to note that 3D-simulations of MHD jet formation have indeed proven 
the feasibility of MHD self-coUimation under the influence of non-axisymmetric perturbations 
in the jet launching region (Ouyed et al. 2003; Kigure & Shibata 2005). Kelvin-Helmholtz 
modes are usually fastest growing and are particularly dangerous for the jet in the super- 
Alfvenic regime. Essentially, it is the "backbone" of the jet flow, i.e. the axial region of 
high field strength and low density which stabilizes the jet against the instabilities as this 
highly sub-Alfvenic region of (Ouyed et al. 2003; Pudritz et al. 2006b). In general, the 3D 
simulations show that the jets starts as a stable outflow, then destabihzes after about 100 
inner disk rotations before it stabihzes again at the time of 200 inner disk rotations. The 
detailed behavior of the different instability modes depend of course on the flow parameters 
and the general setup of e.g. the non-axisymmetric boundary conditions (see Ouyed et al. 
(2003) for details). 

As an interesting fact, we note that all studies of MHD jet formation so far (time- 
dependent or stationary) which lead to a coUimated jet, the disk magnetic fleld proflle 
satisfies the Spruit et al. criterion /x < 1.3. Examples are the stationary (Blandford & Payne 
1982) solution of a self-similar, cold jet with /j = 5/4, simulations of a collimating jet by 
Ouyed & Pudritz (1997a) applying = 1, or by Fendt & Elstner (2000) treating an initially 
dipolar field distribution = 3. 
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3. Model setup 

We perform axisymmetric MHD simulations of jet formation for a set of different bound- 
ary conditions for the accretion disk magnetic field. We apply the ZEUS-3D MHD code 
modified for magnetic diffusivity (Fendt & Cemeljic 2002). As initial state we prescribe a 
force-free magnetic field within a gas distribution in hydrostatic equilibrium. This is essential 
in order to avoid artificial relaxation processes caused by a non-equilibrium initial condition. 
The gas is "cold" and supported by additional turbulent Alfvenic pressure. 

A Keplerian disk is taken as a boundary condition for the mass infiow from the disk 
surface into the corona and the magnetic fiux. The poloidal magnetic field lines are anchored 
in the rotating disk. A gap extends between the central body and the inner disk radius ri 
implying a vanishing mass load for the jet from this region. Compared to our previous 
studies (Fendt & Elstner 2000; Fendt & Cemeljic 2002), several new features are included in 
the present approach. 

- Our major goal is to investigate a broad parameter range of different accretion disk 
magnetic fiux profiles and disk wind density (i.e. mass fiux) profiles. 

- We apply a non-equidistant numerical grid covering a large spatial domain of 150 x 300 
inner disk radii corresponding to e.g. about (6.7 x 13.3) AU assuming an inner disk 
radius of 10 for protostars. 

- We consider a simple parameterization of turbulent magnetic diffusivity, variable in 
space and time, in agreement with also considering turbulent Alfvenic pressure. How- 
ever, compared to our previous studies the magnitude of magnetic diffusivity is low, 
thus not affecting the degree of jet coUimation. 

In the following we give some details for the new features of our the model setup. 



3.1. Resistive MHD equations 



We numerically model the time-dependent evolution of jet formation considering the 
following set of resistive MHD equations. 



^ + V-(p^) = 0, V-B = 0, — J=VxB, 
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+ V(p + pa) + pV* - = 0, 
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with the usual notation for the variables (see Fendt & Elstner (2000)). 

We apply a polytropic equation of state for the gas, p = Kp^ with the polytropic index 
7- Thus, we do not solve the energy equation (4) and the internal energy of the gas is 
reduced ^ to e = p/(7 — 1). The Ohmic heating term in energy equation (4) has shown to 
be generally negligible for the dynamics of the system (e.g. Miller & Stone (1997)). This 
is usually due to the low levels of dissipation presumed in astrophysical cases. Comparison 
of the compression and dissipation terms shows that this holds also for our study (plasma 
beta /9i ^ 1 and magnetic diffusivity 7] ~ 0.01, see below). On the other hand, simulations 
treating violent reconnection processes in stellar magnetospheres and/or current sheets have 
shown that Ohmic heating might establish regions of hot plasma (e.g. Hayashi et al. (1996)). 
As we are primarily interested in the flow dynamics we neglect Ohmic heating. 

Applying a polytropic equation of state simplifies the numerical challenge considerably. 
We do not believe that this simplification affects the general results of our simulations. 
Clearly, certain number values numerically derived in our study may well change (as e.g. the 
exact degree of coUimation) and their interpretation must be strictly limited to a comparison 
within our parameter range applied. Future studies on this topic will have to consider and 
to compare to a non-polytropic equation of state as well. This is, however, beyond the scope 
of the present paper. 

In addition to the hydrostatic gas pressure p we consider Alfvenic turbulent pressure 
assuming pA = p/Pt with Pt — const. This approach has been applied in a series of papers 
by various groups (Ouyed & Pudritz (1997a), Ouyed & Pudritz (1997b), Fendt & Elstner 
(2000), Fendt & Cemeljic (2002), Vitorino et al. (2003), Ouyed et al. (2003), Pudritz et al. 
(2006a)). As we intend to complement these studies, we have followed once again the same 
approach. In general, turbulent Alfvenic pressure may explain the presence of a cold corona 
above protostellar accretion disks as observationally suggested. We have previously argued 
that the existence of such an additional pressure component is inevitable as turbulent Alfven 
waves will be launched naturally in the highly turbulent accretion disk and then propagate 
into the disk wind and outflow (Fendt & Cemeljic 2002). 



^Note, however, that ongoing work by Clarke and collaborators seems to indicate that relaxation of the 

polytropy assumption may afFc'ct the dynamical evolution in certain domains of the jet, in particular regions 
with shocks or contact discontinuities (Ramsey & Clarke 2004) 
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For the poly tropic index we apply 7 = 5/3. This is a value generally assumed for non- 
relativistic astrophysical polytropic gas dynamics. However, note that with such 7, Ohmic 
heating (see above) and also turbulent Alfvenic pressure is potentially neglected for the 
thermodynamics. Detailed (linear) analysis of magnetohydrodynamic turbulence in media 
under various conditions has shown that sub-Alfvenic pressure exerted by Alfven waves 
generally obeys a polytropic equation of state with polytropic index 7t = 1/2. This has been 
derived analytically for well defined model assumptions^ which are applicable for a variety of 
physical states (see McKee & Zweibel (1995)) and has also proven numerically (Gammie & 
Ostriker 1996; Passot & Vazquez-Semadeni 2003). However, the exact number value of the 
polytropic index for turbulent Alfvenic pressure in the case of expanding and accelerating 
magnetohydrodynamic outflows is not known (if it exists at all). Therefore, for simplicity, we 
apply a turbulent Alfvenic pressure polytropic index equal to that of the gas, 7t = 7 = 5/3. 
The effective polytropic index for multi-pressure polytropes {N pressure components p„) 
would be given by 

[Pn\ A7 + 7t 

(see Curry & McKee (2000) and references therein) indicating that for low /3t the polytrope 
would be dominated by the turbulent Alfvenic pressure term. For 7t = 7 = 5/3 and 
/3t — 0.03, also 7efr ~ 5/3. 

We solve the MHD equations applying the ZEUS-3D code (Stone & Norman 1992a,b; 
Hawley & Stone 1995) in the axisymmetry option for cyhndrical coordinates (r, 0, 2;). We 
added physical magnetic resistivity to the original ZEUS-3D ideal MHD code (see description 
and extensive tests in Fendt & Cemeljic (2002)). Wc normalize all variables to their value 
measured at the inner disk radius r-^. For example, time is measured in units of the Keplerian 
period at the inner disk radius, and for the density p — > p/ pi- The point mass for the 



gravitational potential $ = — l/vr^-F^ is located at the origin. Finally, the normalized 
equation of motion treated numerically is 

The coefficients /3i = ^npi/Bj and 5i = Piv\Jpi with the Keplerian speed VK,i = {GM/riY/'^, 
correspond to the plasma beta and the Mach number of the gas at the inner disk radius 
(Ouyed & Pudritz 1997a) . For a "cold" corona with > we have A = l/((5i(7 - 1)/7 - 1) ■ 



thermodynamic system undergoing temporal changes the adiabatic index of (compressible) Alfven 
waves is instead 7wawc = 3/2 (McKcc & Zweibel 1995). Pwave ^ P^^"^ applies in steady state media or Alfven 
waves propagating in a density gradient. 
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In the following primes are omitted and only normalized variables are used if not explicitly 
declared otherwise. 



3.2. Turbulent magnetic diffusivity 

Magnetic diffusivity caused by turbulent motion in the jet material affects both coUima- 
tion and acceleration of the jet (Fendt & Cemeljic 2002). Magnetic diffusivity also plays an 
essential role in the jet launching process in the accretion disk (Li 1995; Ferreira 1997). The 
magnetic diffusivity rj we apply in Eq. (3) is "anomalous" and is considered to be provided 
by macroscopic MHD instabilities. As the jet launching accretion disk is highly turbulent 
itself, we naturally expect that the turbulent pattern will propagate into the jet when the 
disk material is lifted up from the disk surface. The turbulence evolution along the jet will 
certainly be affected by advection and stretching and may vary in different regimes. How- 
ever, the exact dynamical evolution of turbulence in jets has not yet been investigated and 
is therefore poorly understood. 

Supposing that the magnetic diffusivity is primarily due by the same turbulent Alfvenic 
waves which are responsible for the turbulent Alfvenic pressure applied in the simulations, 
we have derived a toy parameterization relating both effects (Fendt & Cemeljic 2002). Such 
an approach extends and generahzes the studies by Ouyed & Pudritz (1997a), Fendt & 
Elstner (2000), and Pudritz et al. (2006a). We achieve this by parameterizing the turbulent 
magnetic dijfusivity similar to the Shakura-Sunyaev parameterization of turbulent viscosity, 
Vt = o^mvl, where < 1 and I and v are the characteristic dynamical length scale and 
velocity of the system, respectively. 

The turbulent magnetic diffusivity rjt can be related to the Alfvenic turbulent pressure 
Pa by first selecting the turbulent velocity field Vt (i.e. the turbulent Alfvenic velocity) as 
characteristic velocity. This gives rjt — a^aVtl and with Pi — {cs/vtY and Cg = jp/p, if 
follows 

2 7 P ,2 7 /VN 

^^t = 75--, or v^ = —p^ . (7) 

For the chosen polytropic index this implies a normalized magnetic diffusivity rj ~ p^^^ if 
I is constant. The power index for r]{p) would be affected by a change in the polytropic 
index for the turbulent Alfvenic pressure (see discussion above). The characteristic length 
scale / of turbulent magnetic diffusivity in the jet is not a known quantity. If we assume 
the turbulence being generated within the disk and launched from there into the outflow, we 
may assume an initial length scale / corresponding to the Shakura-Sunyaev parameterization 
of the disk turbulent viscosity, i.e. limited by the disk height / < h{r). In our normalization. 
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this would imply I ~ 0.1 for regions close to the star as typically h/r ~ 0.1 and I ~ 1.0 at 

r ~ 10. For regions at greater distance to the star, / could be correspondingly larger. A fully 
self-consistent treatment of Alfvenic turbulence would have to consider also its temporal and 
spatial expansion as well as damping or decay as the turbulence cells propagate along the 
collimating jet. This is clearly beyond the scope of our purpose and we instead assume a 
constant characteristic length scale /. 

In our simulations we apply sub- Alfvenic turbulence, j3t = 0.03 and Si = 100 (Ouyed & 
Pudritz 1997a; Fendt & Ccmcljic 2002) resulting in a magnetic diffusivity r] = rjop^^^ with 
TjQ = 0.03. This would correspond to a number value (ctmO — 0.04. 

We finally note that the magnitiidc of normalized magnetic diffusivity is in the same 
range as typically applied in the literature (Hayashi et al. 1996; Miller & Stone 1997; von 
Rekowski & Brandenburg 2004). In a previous study investigating the relation between 
magnitude of magnetic diffusivity and degree of jet coUimation we have shown that such 
(low) level of diffusivity does not affect the degree of jet coUimation - jet de-coUimation as 
due to magnetic diffusivity applies above a critical diffusivity about a factor of 100 larger 
as used in the present paper (Fendt & Cemeljic 2002). In fact, as our goal is to study 
the relation between jet coUimation and disk magnetic field distribution, the magnitude of 
diffusivity was chosen such that is will not affect the degree of coUimation. 

3.3. Numerical grid 

We use the "scaled grid" option by ZEUS with the grid element size decreasing inwards 
by a factor of 0.99. The numerical grid size is (256 x 256) resulting in a a spatial resolution 
from (0.13 x 0.26) close to the origin to (1.6 x 3.2) at the opposite corner. The physical 
grid size for all simulations is {r x z) = (150 x 300)ri corresponding to (6.7 x 13.3) AU for 
r; ~ 10 R0 if applied for a protostar. This is several times larger than in previous studies 
(Ouyed & Pudritz 1997a; Fendt & Ccmcljic 2002; Pudritz et al. 2006a) and comparable to 
the observational resolution for young stars (Bacciotti et al. 2002). 

3.4. Initial conditions 

Prescribing a stable initial condition is essential for any time-dependent simulation. 
Otherwise, the early evolution of the system will be dominated by artificial relaxation pro- 
cesses of the instable initial setup. This is particularly important if only few evolutionary 
time steps are computed. 
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We start our simulations from a density stratification in hydrostatic equifibrium with an 
initially force-free (and also current-free) magnetic field. Such a configuration would remain 
in its initial state if not disturbed by the given boundary conditions. The initial density 
distribution is p(r, z) = (r^ + ^2)-3/4 (Quyed & Pudritz 1997a; Fendt & Cemeljic 2002). 

The initial magnetic field distribution is calculated using a stationary, axisymmetric 
finite element code described elsewhere (Fendt ct al. 1995). Essentially, this code calcu- 
lates the axisymmetric force-free magnetic flux distribution for any given boundary value 
problem. We compute the 0-component of the vector potential"^ using a numerical grid 
with twice the rcsohition of the grid applied in the ZEUS code. Having obtained the vector 
potential of the initial field distribution by the finite element code, we derive the magnetic 
field distribution for the time-dependent simulation with respect to the ZEUS-3D staggered 
mesh (see Fendt & Elstner (2000). Our approach guarantees |V • -B| < 10~^^, respectively 
iJx 5| = 0.01|S|. 

For the finite element grid, we prescribe the disk magnetic flux ^{r,z = 0) — r'^* as 
Dirichlet boundary condition, corresponding to a poloidal magnetic fleld Bp{r, 0) ~r~^ with 
fjL — 2 — //^. Along the axial boundary ^(r = 0,z) — 0. Along the outer boundaries we 
prescribe a magnetic flux decreasing towards the axes (here the exact flux proflle is marginally 
important as the initial fleld will be immediately perturbed by the outflow). 

3.5. Boundary conditions 

Along the r-boundary we distinguish between the gap region extending from r = to 
r = Tj = 1.0 and the disk region from r = 1.0 to r = Tout- The disk region governs the 
mass inflow from the disk surface into corona (denoted also as " disk wind" in the following) . 
The poloidal magnetic fleld along the r-boundary is fixed in time and is, hence, determined 
by the choice of the initial magnetic fiux distribution. Thus, the magnetic fiux across the 
equatorial plane is conserved, ^(r, 2; = 0;t) = \E'(r, z = 0; t = 0). The poloidal magnetic 
field profile along the disk is chosen as a power law, Bp{r, 0) r^r~^. Equivalent to Ouyed & 
Pudritz (1997a) the disk toroidal magnetic field component is force free with = /li/r for 
r > n, Hi = B^^i/Bi. 



The (^component of the vector potential rA^ is identical to the magnetic flux \t'(r, z) = (l/27r) / Bp-dA 
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3.5.1. Hydrodynamic boundary conditions 

The hydrodynamic boundary conditions are "inflow" along the r-axis, "reflecting" along 
the symmetry axis and "outflow" along the outer boundaries. Matter is "injected" from the 
disk into the corona parallel to the poloidal magnetic fleld with low velocity ■Uinj(r, 0) = 
J^iVK{'>^)Bp/Bp and with a density pinj(r, 0) = pip(r, 0). In some simulations we applied a 
slightly enhanced mass flow rate from the outermost grid elements of the disk by increasing 
the inflow velocity by a certain factor. This choice helped to deal with some problems at the 
intersection of inflow and outflow boundary conditions just in the corner at (rmax, = 0), 
and does not affect the flow on larger scales. Along the gap, the mass flux into the corona 
is set to zero. We chose a power law for the inflow density proflle, 

Pi„j(r,0) = po(r)=Pir-^''. (8) 

As a flrst choice, it is natural to assume /ip ~ 3/2 as this is the disk density proflle both 
for the Shakura-Sunyaev accretion disk and for advection dominated disks (see e.g. the 
self-similar Blandford & Payne solution, or the Ouyed & Pudritz (1997a) model). This is, 
however, the average density profile of a thin disk and one may expect deviations from that 
profile for several reasons. Those may arise (i) as the structure of jet-driving disks can be 
altered by the outflow, (ii) due to magnetic flelds close to equipartition, or may be simply 
suggested by the fact that (iii) we do not yet fully understand the mass loading of the jet / 
disk wind out of the accreting matter. We have therefore applied different density proflles 
lip = 0.3, 2.0 for the disk wind mass flux in our simulations. 

For the injection velocity, one may assume that the initial disk wind speed is in the 
range of the sound speed in the disk, Vi^^j(r) ~ Cs(r) ~ Wkep ~ r~-^/^. However, the actual 
velocity along the disk surface will be different from the thin disk estimate and will, moreover, 
depend in detail from the launching mechanism. For simplicity, all simulations presented in 
the paper assume the same injection velocity. Therefore, the total mass flux profile is varied 
only by the density profile. 

3.5.2. The disk magnetic field profile 

The disk magnetic flux profile is less constrained, or better, less understood. A thin disk 
equipartition magnetic fleld will follow a proflle ~ r~^^* (Blandford & Payne 1982). How- 
ever, it is not clear to date where the disk/ jet magnetic field originates, how it is generated 
and how it evolves. We find it therefore valuable to explore a wider range of disk magnetic 
fiux profiles. This is the main goal of the present paper. The choice of having a certain disk 
magnetic field distribution fixed in time should be understood as this field structure being 
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determined by the internal disk physics (implying that different disk internal conditions may 
lead to a different disk flux distribution). In this respect, we do not care about how the 
magnetic field has been built up in time. Instead, we assume that internal disk processes like 
diffusion, advection, or turbulence work together in generating and supporting a large scale 
stationary field structure with the chosen disk surface field distribution. In the introduction 
we have already stressed the point that concerning the internal disk physics and the inter- 
action between disk and outflow (in both directions) some of the parameter runs considered 
may finally be proven to be more realistic than others if compared to a fully self-consistent 
simulation including the disk internal disk physics and dynamics in relation with the out- 
flow itself. So far, this problem has not been solved in generality, although there are first 
simulations treating the disk-jet system over many hundreds of rotational periods applying 
simplified disk models (Casse & Keppens 2002; Kudoh et al. 2002; Ccmeljic & Fcndt 2004; 
von Rekowski & Brandenburg 2004). In particular we mention the work by (von Rekowski 
& Brandenburg 2004) who in their simulations did not constrain the disk magnetic field but 
let it evolve by a dynamo process from an initial stellar magnetosphere. Indeed the disk 
magnetosphere evolves into an open field structure, however the radial profile of the disk 
field is not explicitly mentioned. The outfiows from disk and star are pressure driven but 
a magneto-centrifugally driven outfiow region co-exists originating in the inner part of the 
disk. However, no collimated jet fiows have been observed in these simulations covering a 
area of 0.1 x 0.2 AU. 

Some insight into the parameter space investigated in our simulations may be gained 
by comparing them to stationary models of accretion-ejection structure. As an example, in 
the following wc discuss the solutions presented by Ferreira (1997). Ferreira et al. apply 
self-similar and stationary MHD for the disk-outflow system and derive the so-called ejection 
(or mass loading) index ^ = d\n Mg^cc/dlnr as the leading parameter governing the structure 
of the outflow. They flnd this parameter being constrained by a lower limit of 0.004 due to 
the vertical disk equilibrium, and a maximum value ^ < 0.15 due to the condition of having a 
trans-Alfvenic jet^. Since by definition ,^ > 0, it follows that Hp < 3/2. In fact, for fXp = —3/2 
the Ferreira et al. approach gives the Shakura-Sunyaev disk density distribution and ^ = 0, 
i.e. a constant accretion rate dMs^{r) and no outfiow. Their upper limit ^ = 0.15 would 
imply a very narrow range for the radial density profile, 1.35 < /ip < 1.5. However, in our 
simulations wc obtain trans-Alfvenic outflows also for lower Hp. Our interpretation is that 
this constraint might be an artifact of the self-similar approach undertaken by Ferreira et al. 
Without considering the disk structure self-consistently in a time-dependent simulation this 
question cannot be answered. 



^Note that with dMa^c/dr ~ dMjet/dr we have iip = 3/2 — ^ 
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The self-similar study further provides a constraint for the profile of the magnetic flux 
configuration < /x^ < 2 (note that j3 in Ferreira's notation corresponds to j.!^ in our's). 
This is in agreement with all parameter runs in our study. Another relation found by Ferreira 
et al. is that ^q, can be related to ^hy p^, = 3/4 + ,^/2, or = 1/2 + /ip/2, respectively, 
(the number value 3/4 indicates the Blandford-Payne solution). This quantifies the relation 
between the profile of the disk mass loss and that of the disk magnetic field. The above 
mentioned constraint /Xp < 3/2 implies that // < 5/4, which tells us that our simulations in 
the upper part of Tab. 1 are inconsistent with the self-similar disk- jets discussed by Ferreira 
et al. However, considering the fact that (i) the steepness of the field profile is not so different 
from the other runs and that (ii) the analytical constraint is limited by the self-similarity 
assumption, we still think that these runs provide valuable scientific information. 

Finally, we can only stress the point once more that both the the mass loading from 
disk into the outflow and the field structure will eventually defined by the disk physics and 
the constraints to it set by the existence of an outflow. 



3.6. Jet flow magnetization 

In ideal MHD the magnetic field is "frozen" into the matter which itself slides along 
(but not parallel to) the field lines. In an axisymmetric and stationary configuration the 
magnetic field lines are located on magnetic flux surfaces ^(r, z) and the velocity vector can 

— * 

be expressed SiS v — k{^)B/ p — rVt-pi^)-, where k{^) = {dM /d^) — {pvp/ Bp) is the mass 
flow rate per magnetic flux surface and Of(^) the iso-rotation parameter^. One of the most 
important parameters for MHD wind theory is the magnetization parameter Michel (1969), 

a = P. ^ , (9) 

where M(^') is the mass flow rate within the magnetic flux surface The magnetization 
is inverse to the mass load a ~ k~^. For outflows expanding in spherically radial direction, 
Michel (1969) derived an analytical relation between the asymptotic outflow velocity and 
the magnetization, v^o ~ cr^^^ (Michel scaling). The outflow coUimation results in a varia- 
tion in the power index of the Michel scahng (Fendt & Camenzind 1996). Essentially, the 
magnetization summarizes the important parameters for launching a high velocity outflow - 
rapid rotation, strong magnetic fleld and/or comparatively low mass load. It will therefore 



^For simplicity, the iso-rotation parameter flp can be interpreted as the angular velocity of a magnetic 
field line. Both fc(\l/) and f2F(^') are conserved quantities along the magnetic flux surface. 
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be interesting to study also the influence of the magnetization profile across the disk wind on 
the coUimation of this wind into a jet. Having prescribed a power law for the profiles of the 
disk wind initial poloidal velocity Viajir) = t'inj'yKep(^), the rotation ^^{r) = l^Kep ~ r~^/^, 
the density and magnetic field profile, the disk wind magnetization profile cro(r) = cr(r, ^ = 0) 
follows a power law as well, 

p{r) 

(hence, = l^p — '^1^ — 1/2). A choice of parameters as Piniir) ~ r~^/^, Bp{r, z = 0) ~ 
(Ouyed & Pudritz 1997a; Fendt & Cemeljic 2002) will result in a radial disk wind 
magnetization profile (7o(r) ~ r~^. The self-similar Blandford & Payne model gives pa — 
—3/2. In a recent paper Pudritz et al. (2006a) have followed an approach similar to the 
present paper, varying the mass load k{r) at the disk surface as k{r) ~ r~^, ^-3/4:^ r~^^'^, r~^l^. 
This would correspond to magnetization profiles as \Xfj — 1,0.75,0.5,0.25, respectively, a 
parameter space which is also covered in our paper. 



4. Results and discussion 

We have run a large number of MHD jet formation simulations covering a wide parameter 
range concerning both the disk wind poloidal magnetic field profile and the "injection" 
density profile along the disk surface (see Tab. 1). The results presented here are preliminary 
in the sense that not all simulation runs could be performed over time scales sufficiently long 
enough for the MHD outflow as a whole to reach the grid boundaries or to establish a 
stationary state. This is due to the large grid size in combination with steep gradients 
in the disk wind parameters which in general allow only for a weak outflow from large disk 
radii. Nevertheless, our results allow for a flrm interpretation of the overall outflow structure 
and evolution in the MHD jet formation region. In most cases we have obtained a quasi- 
stationary state for the outflow within the inner region (over at least 50% of the grid). In 
some cases the nature of the outflow prevents that a stationary state can be reached. 
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Fig. 1. — Simulation run ill. Mass density in grey colors (logarithmic scale between between 
the density values shown at the color bar on the very top) and poloidal magnetic field lines 
{black lines) at times t = 0, 100, 500, 1500, 2090 (from top to bottom). The r and 2;-axis point 
in vertical and horizontal direction, respectively. The dynamic range of the color coding 
is increased by setting densities above (below) the values given at the color bar to black 
(white) . The thick dark grey and light grey contours indicate the run of the Alfven surface 
and the fast magnetosonic surface, respectively. 
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Fig. 2. — Simulation runs with different disk boundary conditions (see Tab. 1). Density 
and poloidal magnetic field lines for p4 {t = 2370), il5 {t = 5000). p20 {t = 1810), a2 
{t — 2000), ilO {t — 3000) from top to bottom. See Fig. 1 for further notation. 
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Fig. 3. — Simulation runs with different disk boundary conditions (see Tab. 1). Density and 
poloidal magnetic field lines for pl4 {t = 1840), i2 {t = 4000), i3 {t = 4000), pl6 (t = 2790), 
clO {t — 1400) from top to bottom. See Fig. 1 for further notation. 
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4.1. The general flow evolution 

The main features of jet formation in the model scenario apphed have been described 
in the previous hterature (Ouyed & Pudritz 1997a, 1999; Fendt & Cemeljic 2002) and will 
not repeated here in detail. In summary, the early evolution (t < 100) is characterized by 
the propagation of torsional Alfven waves launched by the differential rotation between disk 
and corona (see Fig. 1). The Alfven waves shghtly distort the initial field structure into a 
new state of hydrostatic equilibrium which will remain stationary until reached by the MHD 
outflow launched from the disk. The region beyond the wave front remains undistiirbed. 
Torsional Alfven waves are still propagating across the grid, when the MHD wind begins 
to accelerate from the disk surface and becomes collimated beyond the Alfven surface. The 
outflow rams into surrounding magnetohydiostatic corona developing a bow shock which 
propagates with velocities comparable to the material speed. The flow structure which is 
left when the bow shock has passed the computational grid is a pure disk wind/jet, generated 
and energetically maintained by the disk rotation, the continuous mass injection from the 
disk and the disk magnetic flux. We note that at these later evolutionary stages the outflow 
has completely swept out of the grid the initial condition and is governed solely by the 
boundary conditions. 

As a general result, outflows launched from a flat disk magnetic field profile (// < 0.5) 
do not reach a steady state. Instead, a wavy pattern evolves close to the symmetry axis 
with almost constant mass flow rate in axial direction but variable mass flow rate in lateral 
direction. This case is similar to the extreme example of an initially vertical magnetic held 
examined by Ouyed & Pudritz (1997b). 
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Table 1: Summary table of our simulation runs. Shown are parameters for the disk poloidal 
magnetic field profile fi, the disk wind density profile fip, and the disk wind magnetization 
profile /icT- The inner disk radius poloidal field strength is Bp^i. The density at this radius is 
p^j = 1, but is sometimes lowered by a factor of two (denoted by a * in column 1). Mass flow 
rates in r- and 2;-direction along the three sub-grids {rmax,i x -^max.i) = (20.0 x 60.0), (40.0 x 
80.0), (60.0 X 150.0) are denoted as Mz,i, Mr,h respectively. The average degree of coUimation 
measured by the mass flow rates in z and r-direction, and normalized by the area threaded, 
is <C>- The maximum flow velocity obtained outside the axial spine is Umx- 
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Fig. 4. — Degree of collimation measure (denoted by <C>) as defined in Sect. 4.2 and power 
law index of the disk poloidal magnetic field profile fi (see Tab. 1). The three sub-figures 
refer to simulations applying inflow density profiles, /Xp = —1.5 (top), /ip = —1.0 (middle), 
= —0.5 (bottom). Data points with bars indicate a variable collimation degree resulting 
in most cases from axial instabilities of highly coUimated outflows. 
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4.2. How to measure the degree of collimation? 

The main goal of this study is to investigate the relation between the accretion disk 
magnetic flux proflle and the collimation of the outflow. In order to quantify such an rela- 
tion we need to define a measure for the degree of outflow collimation. As a matter of fact, 
some ambiguity exists about such a definition and a few words of clarification are needed. 
What is usually observed as jet fiow, is an elongated feature of radiation emitted in combi- 
nation with a signature of high velocity. This is then associated with (i) an elongated mass 
distribution ("collimation") emitting the radiation and (ii) a strong axial velocity component 
("acceleration"). Therefore, from the theoretical point of view, it is natural to apply the 
directed mass flux as measure of jet collimation (see below). 

In this respect, we also note the different meaning of the "degree of collimation" and 
the "ability to coUimate". The first term refers to the status of the established - coUimated 
- outfiow as directed mass flux which is probably equivalent to the observed jet. The second 
term refers to the question of how well an initially weakly collimated disk wind can be turned 
into a collimated stream by means of internal MHD forces. It is straight forward to expect 
that for the same disk magnetic field profile, a more concentrated disk wind density profile 
(i.e. large /Xp) will result in a more concentrated density profile of the asymptotic jet flow 
as well (presumably leading to a smaller jet radius observed). However, both a wide and a 
narrow jet, i.e. a jet with a flat or steep asymptotic density distribution, may have the same 
asymptotic opening angle (e.g. a cylindrical shape). 

For ideal MHD jets in stationary state, collimation can be measured by the opening 
angle of the magnetic field fines. In the case of resistive jets, or if the outflow has not yet 
reached a steady state, this is not appropriate, as mass flux direction and magnetic fleld 
direction are not aligned. 

In the present paper, we consider the directed mass flux as appropriate quantitative 
measure of outflow collimation. For simplicity, we measure the degree of collimation by the 
ratio C of mass flux in axial versus lateral direction, 

C^^^ 2vr/;-rpMr ^^^^ 
M, 27rrmax Jo pVrdz. 

This parameter depends on the size of the integrated volume. We determine the mass flow 
rates applying three differently sized volumes, Mz,i, -^r,i, considering different sub-grids of 
the whole computational domain of 150 x 300 n corresponding to cylinders of radius and 
height (r^ax x z^ax) = (20.0x60.0), (40.0x80.0), (60.0x150.0) (see Tab. 1, Tab. 3, Tab. 2). 
We then consider the mass flux ratio ( normahzed to the size of the areas passed through by 
the mass flow in r- and z-direction, d, ^2, Cs- The ratio of grid size (rmaxZ-^max) as displayed 
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in our figures converts into a ratio of cylinder surface areas of A^/A,- ~ 0.5(rinax/-Zmax)- Thus, 
a mass flux ratio C2 = 2 as defined by Eq. (11) corresponds to a mass flux ratio normalized 
to the surface area threaded of ^2 = 8 which could be considered as "very good" coUimation. 

The difference in degree of coUimation calculated for differently sized volumes, Ci, C2, Cs) 
may also indicate different evolutionary states in different regions of the outflow. At a certain 
time, the oTitflow may already be well coUimated close to its origin and also have reached a 
stationary state in this region. However, on larger spatial scales the same outflow may still 
interact violently with the ambient initial corona (see e.g. simulation runs c3, p4). Some 
simulations could be performed until an almost stationary flow on the whole numerical grid 
has developed (e.g. c8, p8). In other cases a quasi-stationary state of the outflow within the 
inner region (say 50% of the grid) has been reached. In most cases the initial corona has 
completely swept out of the grid. Prom the mass flow rates on the sub-grid we derive an 
average degree of coUimation <C> considered to be typical for that simulation, but depending 
on the evolutionary state of the flow. 

We note another interesting consequence of the large grid applied compared to previous 
studies on smaller grid size (Ouyed & Pudritz 1997a; Pudritz et al. 2006a). In some of 
our simulations which reaches later evolutionary stages of the outflow the Alfven surface 
leaves the grid in axial direction (sec e.g. run 13). Thus, these outflows are sub-Alfvenic 
in their outer layers and the jet structure will consist of three nested layers (cylinders) of 
different magnetosonic properties - an outer sub-Alfvcnic layer, a super-Alfvenic / sub-fast 
magnctosonic layer at intermediate radii, and a super-fast magnetosonic jet in the inner 
part (here we neglect the very inner axial spine where there is no mass inflow from the 
disk but potentially from a stellar wind). The extension of these layers depends on the 
parameters /i and /ip and will be discussed in more detail below. Por example, for simulation 
run ilO, we estimate the asymptotic radius of the Alfven surface to about 120ri, and the 
asymptotic radius of the fast magnetosonic surface to about 50ri. Por comparison, in the 
original approach (Ouyed & Pudritz 1997a) passes the radial grid boundary at a height of 
about z = 60. A multi-layer structure of jets in respect to the magnetosonic surfaces has 
been proposed in stationary studies by Fendt & Camenzind (1996), however, then suggesting 
layering inverted to what we flnd now (ie. a super fast-magnetosonic flow in the outer part 
respectively) . We expect that such a stratiflcation in lateral direction should be essential for 
the jet stability. This particular issue has not yet been treated in the literature. 
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Fig. 5. — Degree of coUimation (denoted < C >) ^ defined in Sect. 4.2 and power law- 
index of the disk wind magnetization profile (see Tab. 1). This figure combines all three 
sub-figures of Fig. 4 applying /i^ — /ip — 2/1 — 1/2. Data points with bars indicate a time 
variable coUimation degree resulting in most cases from axial instabilities of highly coUimated 
outflows. 
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4.3. Degree of coUimation - MHD simulations 

Simulation run ilO with fi = 1.0, fip = 1.5 serves as our "reference run" resembling 
just the same parameter set as applied in the exemplary solution by Ouyed & Pudritz 
(1997a), but now performed on a substantially larger physical grid. After t = 3000 (i.e. 
3000 inner disk orbits) a stationary outflow is obtained over a large fraction of the grid. The 
run of the Alfven surface at large radii indicates that the outermost layers of the outflow 
have not yet reached a stationary state (see Fig. 2). We obtain degrees of coUimation of 
Ci,2,3 = 9.0, 9.2, 12.5 resulting in an "average" degree of coUimation of <C>= 10 (see Tab. 1). 
As the outflow has not yet settled in a stationary state over the whole grid, these values 
will vary slightly for later times. Note that even 3000 orbits of the inner disk correspond to 
about only two orbits at the outer disk radius at 150 Ti. The sequence Ci,2,3 also indicates 
an increasing degree of coUimation along the outflow. 
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Fig. 6. — Simulation runs with different disk boundary conditions (see Tab. 1). Impact of 
a total disk magnetic flux variation (factor 4). Density and poloidal magnetic fleld lines for 
pl8 {t = 2210), pl9 {t = 1640) from top to bottom. See Fig. 1 for further notation. 
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4.3.1. Disk magnetic field profile and collimation 

We now compare simulation runs with different disk magnetic flux proflle fi,fi<!j, but 
similar inflow density proflle fip. We have run about 40 different simulations which in total 
required about one year of CPU time on various work stations. 

Table 1 shows the main input parameters for the each simulation run, in particular the 
power law index of the magnetic fleld and density distribution. Table 1 further shows the 
main parameters of the resulting outflow, namely the mass flow rates in lateral and axial 
direction defining the collimation degree. What is also indicated is the physical time step 
T at which the mass fiow rates are calculated. This is in most cases, but not always, the 
final time step of the simulation. When comparing different simulation runs, one has to keep 
in mind that they evolve with different dynamical speed. Thus, if possible, only outflows 
at late (or equivalent) evolutionary stages will be compared, which may have been reached 
after different physical time (as e.g. stronger magnetized outflows evolve faster, see below). 

In general, we find an increasing degree of collimation with decreasing slope of the 
disk magnetic field profile^ (see Fig. 4). This holds for all of the disk wind density profiles 
investigated so far. A steep density profile generally leads to a higher degree of jet collimation 
which is, however, not surprising as the mass fiux is more concentrated along the axis just 
by definition of the boundary condition. 

For the simulations with relatively flat density proflle {fip ~ 0.3, 1.0), a strong increase 
in collimation degree can be observed just below n = 0.6. However, simulations with a flat 
disk magnetic fleld proflle (below /i = 0.5) do not reach a stationary state within the time 
scale of our simulations. In general, the degree of collimation for these outflows is high. The 
axial mass flow rate stays more or less constant whereas the radial mass flux varies in time, 
resulting in a "wavy" pattern evolving along the outflow axis. This has been discussed by 
Ouyed & Pudritz (1997b) for the case of an initially vertical coronal magnetic field and in 
the context of knot generation in protostellar jets. 

Now we compare our results obtained for the different density profiles quantitatively. We 
find that a physically meaningful classification can be achieved by shifting and spreading the 
three figures (Fig. 4) not in vertical but in horizontal direction. Essentially, this corresponds 
to a comparison collimation degree versus magnetization by applying /i^ — /ip — 2/1 — 1/2. 
Indeed, the resulting diagram (Fig. 5) shows a convincing correlation between the magne- 



^Thc mass flow rates in simulations with flat magnetic field profile n < 0.4 are difficult to determine as 
those outflows do not reach a stationary state and show a time-dependent variation in the lateral mass flux 
close to the axis 
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tization profile power law index ji^r and the degree of coUimation obtained <C>- As we 
naturally expect that the inflow parameters will change during the life time of an accretion 
disk, our results might explain why not all accretion disks do launch a jet and why disks do 
not launch jets for all their life time. 

The above mentioned steep rise in coUimation degree for < 0.6 is now indicated for a 
magnetization profile ix^ > —1.0. Obviously, the slope of this rise depends in particular on 
the reliability of simulation runs a6, c8 and i2. We note that simulation i2 and c8 follow a 
relatively fiat density profile, and did evolve up to several 1000 disk orbital periods into an 
almost stationary state on a global scale. Simulation a6 remains sub-Alfvenic over most of 
the grid size and has developed an instable flow pattern along the jet axis. 

The width of the {fi^^ - <(^>)-correlation shown in Fig. 5 is due to detailed differences 
in the parameter setup. The data points shown for a certain simulation which is deviating 
from the general trend of the correlation indeed have a particular origin. For example, 
simulation run c3 is strongly magnetized with the Alfven surface very close to the r-axis. 
The same holds for simulation run 14. Therefore, both simulations do not consider a typical 
Blandford-Payne magneto-centrifugally launched jet, but have the Alfven points close to the 
foot points of the field lines (hence a short lever arm). However, in summary we conclude 
that the numerical data give a clear correlation between the power law index of the disk 
wind magnetization profile and the degree of coUimation. 
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Table 2: Summary of simulation runs obeying an artificial decay of the toroidal field compo- 
nent. The parameter Jbs gives the factor by which the toroidal magnetic field is artificially 
decreased each time step (see Tab. 1 for other notations). 



run 






Mo- 




T 


/S3 


Mz,l, Mr,l 




Mz,a, Mr,a 


Cl,C2,C3 




<c> 


plO 


1.5 


0.8 


-2.7 


0.40 


1800 


0.1 


-.04, 1.58 


0.07, 2.45 


0.11, 3.19 


(0.2), 0.1, 0.2 


0.18 





pll 


1.5 


0.8 


-2.7 


1.20 


3100 


0.1 


-.051,1.59 


0.057, 2.47 


0.064, 3.40 


(0.2), 0.1, 0.1 


0.17 





pl7 


1.0 


1.5 


-1.0 


0.232 


1890 


0.9 


0.30,1.08 


0.82,0.97 


0.96,1.39 


1.7, 3.4, 3.4 


0.17 


3 


il2 


1.0 


1.0 


-1.5 


0.387 


1520 


0.9 


0.60,2.65 


1.73,3.20 


2.59,5.66 


1.4, 2.2, 2.3 


0.29 


2 


il3 


1.0 


1.0 


-1.5 


0.246 


1020 


0.9 


.092,3.19 


1.33,3.73 


3.48,7.38 


0.17, 1.4, 2.4 


0.32 


2 


pl2 


0.8 


1.5 


-().() 


0.922 


870 


o.y 


0.87, O.liO 


1.5.'), 1.19 


1.20, 1.75 


7.0, 5.2, 8.0 


0.27 


1 
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4-3.2. Field strength and collimation 

While the slope of the disk magnetic flux / magnetization proflle clearly governs the 
asymptotic collimation of the outflow, the amount of total disk magnetic flux is seemingly 
less important. We have investigated different absolute magnetic fluxes for the same flux 
profile and did not detect any relation with the degree of collimation. Simulations with high 
magnetic field strength are CPU time consuming due to the Alfven speed time stepping. On 
the other hand, the jet evolves faster and also propagates across the numerical grid within a 
shorter period of time (physical units). We thus reach similar evolutionary stages at earlier 
physical time (but after more numerical time steps). 

For reference we compare simulation runs il4 and il5 with field strength Bp^i{il5) — 
0.5Bp^i{ilA) and initial fiow density pi„j(il5) = 0.5pi„j(il4); runs pl3, pl4 and al with 

Bp^i{pl3) = 3.0Bp^i{pl4:) = 3.73Bp^i{al) and pinj{il5) = 0.5pjnj(^14); and runs pl8, and pl9 
with Bpi{pl9) = 2.0Bpi{pl8) (see Fig. 6). Table 1 shows that the collimation degree is 
indeed similar for each example, as it is the evolving jet structure (see examples pl8 and 
pl9. Fig. 6). On the other hand, the total magnetic flux governs the asymptotic velocity 
gained by the jet (see below). 

4.4. Jet collimation for highly diffusive toroidal field 

As mentioned above (Sect. 2. 3), studies in the literature have proposed a "poloidal col- 
limation" of the jet, i.e. jet collimation due to the pressure of an ambient poloidal magnetic 
fleld (Spruit et al. 1997; Matt et al. 2003). In order to investigate the feasibility of such 
a process, we have performed a set of simulations applying an artiflcially strong decay of 
the toroidal field component, invented to mimic the toroidal field decay due to the kink- 
instability (Spruit et al. 1997). 

Our simulations (see Tab. 2) follow a parameter setup equivalent to those in Tab. 1 
with the single exception that the toroidal field is artificially decreased by either a factor 
/b3 = 0.9 or 0.1 at each numerical time step^. Thus, this part of our paper considers a "toy 
model" mimicking an artificially enhanced magnetic diffusivity, only affecting the toroidal 
field component^. A fully self consistent approach would involve a tensorial description of 



^This is done in the ZEUS subroutine "ct" treating the constrained transport and updating the magnetic 
field components 

^We note that the artificial toroidal field decay sometimes causes problems for the usual outflow boundary 
condition in axial direction. However, so far the internal structure of the outflow has not been not affected 
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magnetic diffusivity with strong diffusivity for those components affecting the toroidal field 
component. This is beyond the scope of the present paper. 

During the long-term evolution of these outflows the resulting toroidal magnetic field 
is decreased by 5 orders of magnitude compared to the simulations presented before (inde- 
pendent of /bs). The simulation runs we have to compare are pl7 with p2 and ilO, runs il2 

and il3 with ill, run pl2 with pl4 and pl3, and runs plO and pll with il6. All outflows 
with decaying toroidal magnetic field are substantially less collimated (Tab. 2). This is not 
surprising as the collimating tension force by the toroidal field has just been switched off. 
The exception is simulation run pl2 which, due to its steep disk wind density profile, shows 
relatively good coUimation comparable to the previous simulations runs without toroidal 
field decay. 
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Fig. 7. — Toroidal velocity profile across the jet ai z = 160 for selected parameter runs at 
the final stages of simulation. i3 {t = 4000), il {t = 3000), p8 {t = 3200), ill {t = 2000) 
from top to bottom. 



4.5. Jet velocity &c Michel scaling 



In general we obtain a narrow range for the maximum velocity of the super-Alfvenic 
flow at large distance from the origin. The maximal jet velocity is typically in the range 
of the Keplerian speed at the inner disk radius, in agreement with the literature (Ouyed 
& Pudritz 1997a; Casse & Keppens 2002; Kudoh et al. 2002). However, in our simulations 
we noticed a slight trend in the jet velocity in relation to the magnetization a. In order to 
further investigate the Michel scaling for coUimated jets, we have performed simulation runs 
with identical initial field and mass flow profile, but different absolute number values for 
mass flow rate and magnetic field strength - either varying the initial magnetic field strength 
or the inflow density by a certain factor. 

These simulations are summarized in Tab. 3. In particular, we may compare simulations 
p21 to 19, c3 to p4, il4 to il5, pl3 to pl4 and al, il to p7 and c5, and pl8 to pl9 (see Fig. 6). 
Table 3 shows the ratio of the maximum poloidal velocity for either of the simulations (labeled 
as 'A' and 'B') as expected from the original Michel scahng 



and each of the numerical values involved (iso-rotation and injection velocity are the 
same for all examples). The jet velocity has been measured along a slice across the jet at 
z = 160. The poloidal velocity profile has a maximum close to the axis. For each sample, 
the velocity measure of the asymptotic flow has been taken at similar spatial location and 
comparable evolutionary time step. This is an important point as we compare simulations 
which different evolutionary time scales where in some cases the large scale flow has not yet 
settled in a steady state. 

The situation for the toroidal velocity is less pronounced. The typical rotation velocity 
at {z = 160, r = 80) is 1/10 of the Keplerian speed at the inner disk for most of the cases 
investigated. Examples for toroidal velocity profiles across the jet are shown in Fig. 7. 

Our large-scale numerical grid allows for direct comparison of the numerically derived 
velocity structure with recent observations indicating rotation in the jet from the young 
stellar object DGTau (Bacciotti et al. 2002). The numerical grid covers about 7 x 14 AU 
which is about the size of the innermost slit position shown in Fig. 7 of Bacciotti et al. (2002). 
These authors measure radial velocities up to 400 km/s and find two major outflow com- 
ponents - a low velocity and a high velocity component. It is the low velocity component 
of de-projected^'^ jet velocity of about 90 km/s which shows radial velocity shifts of up to 



Assuming an inclination of 38°, see (Bacciotti et al. 2002) and references therein 




1/3 



2/3 



1/3 



(12) 
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12km/s. The radial velocity differences across the jet have been widely interpreted as indi- 
cation of rotation although this is not yet completely confirmed. The observed velocity data 
corresponding to our numerical grid (equivalent to "Slit S3/S5" and "Region I" in Bacciotti 
et al. (2002)) are 64km/s in poloidal velocity (in [Sll]) and 8km/s as total radial velocity 
shift. The Keplerian speed at the inner disk radius for DG Tau is 

assuming a central mass of 0.67 Mq (see Bacciotti et al. (2002) and references therein). 
According to such a normalization, our simulations which fit best to the low velocity jet of 
DGTau are those with a fiat magnetic field profile ~ 0.5 and with moderate field strength, 
namely il, p8 (see Tab. 1) These runs also show a high degree of coUimation which is 
observed as well. The numerically derived rotational velocity for these examples is about 
O.Ivk,! ~ lOkm/s which is surprisingly close to the observed value (see Fig. 7). The launching 
radius can then derived just by following back the streamlines along the coUimating jet flow. 
In the runs il, p8 the launching region is distributed over a relatively large disk area of 
r^SO corresponding to about 3 AU. Note that in these simulations only the inner part of the 
jet is super-Alfvenic. Note also that the toroidal velocity in the numerical models decreases 
with radius in apparent contradiction with the observations of DG Tau (which, however, do 
not reolve the region close to the axis) . We like to emphasize that in a self-consistent picture 
also the high velocity component has to be explained by the same simulation. As another 
note of caution we mention that all MHD simulations of jet formation done so far have shown 
a terminal jet speed of the order of the Keplerian speed at the launching region (Ouyed & 
Pudritz 1997a, 1999; Fendt & Cemeljic 2002; Vitorino et al. 2003). It is therefore difficult to 
interpret the very high velocity 450km/s features observed for DGTau in such a framework 
unless we assume very strong magnetic flux respectively a very low mass loading. 

We finally discuss briefly the outflow velocities obtained in simulations with strongly 
diffusive toroidal fleld. As seen from Tab. 2, the typical velocities in such outflows are 
about a factor 5 lower compared to the standard approach. In particular, this follows from 
comparing simulations p2 to pl7, ill to il2 and il3, and pl3 to pl2. This is not surprising 
as we have artificially reduced the accelerating Lorentz force component (along the field), 
-^L,|| ~ j± X B^, where j± is the poloidal electric current density component jp ~ V x 
perpendicular to the magnetic field. 
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Table 3: Subset of Tab. 1, summarizing simulation runs with similar field and density profiles 
but different relative magnetic flux. For the notation of the jet velocity ratio TZy = Voo,a/voo,b 
and the magnetization ratio TZ^ = (Ja/ctb we refer to the text and Eq. 12. 



run 


1^ 




Pi 




<c> 


7^y^ 


7^v 


c3 


1.5 


0.214 


1.0 


1.07 


15 


1.31 


1.63 


p4 


1.5 


0.321 


1.0 


1.74 


10 






il4 


1.25 


0.355 


1.0 


1.53 


10 


1.26 


1.25 


il5 


1.25 


0.177 


0.5 


1.22 


10 






pl3 


0.8 


0.922 


1.0 


2.33 


15 


2.08 


1.98 


pl4 


0.8 


0.307 


1.0 


1.18 


15 


1.17 


1.09 


al 


0.8 


0.247 


1.0 


1.08 


15 






il 


0.5 


0.112 


1.0 


0.7 


30 


1.59 


1.23 


p7 


0.5 


0.225 


1.0 


0.86 


20 






c5 


0.5 


0.225 


1.0 


0.78 


30 






pl8 


0.5 


0.112 


1.0 


0.5 


5 


2.13 


2.46 


pl9 


0.5 


0.225 


1.0 


0.95 


6 






pl6 


0.5 


0.225 


1.0 


1.1 


7 
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4.6. The run of the Alfven surface 

Another parameter which measures the relative strength of the magnetic field is the 
Alfven Mach number Ma, 

Ml ^ (14) 

At the inflow boundary (= disk surface) the Alfven Mach number varies as 

Mloir) ~ r^i'-i'''-' ~ r-"'^-^/'' = r"'*^, (15) 

with Ha = l^a + 3/2. A classical Blandford-Payne type MHD jet is launched as slow sub- 
Alfvenic disk wind, i.e. Ma,o(?^) < 1 for all radii along the disk surface. In turn, this condition 
constrains the disk wind poloidal magnetic field and density profile, respectively, /Xp > 2/1+1. 
For smaller /i, a situation may occur that at a certain radius the disk wind will be launched 
as super- Alfvenic fiow. The location of this radius depends on the actual values of Sp i and 
Pinj. Therefore, in general sub- Alfvenic disk winds require a positive exponent for the MA(r) 
profile. In the self-similar Blandford & Payne approach we have jifj = —3/2, while in the 
Ouyed & Pudritz (1997a) simulations fi„ = —1 A dipolar (stellar) magnetic field distribution 
with /ip ~ 3/2 implies fj,^ = —1 and, thus, Ma,o('") ~ r^^^^, that means a critically negative 
exponent. The disk wind launched will be super- Alfvenic (Fendt & Elstner 2000). 

Note that the relation between the magnetization a and the collimation measure <C> 
shown in Fig. 5 holds similarly for the Alfven Mach number as a*a = A^o- + 3/2. 

Figure 8 summarizes our results concerning the shape of the Alfven surfaces in respect 
of outfiow collimation degree. The showcase plot indicates the run of the Alfven surface of 
different simulations, overlaid by the numerically observed degree of collimation for these 
outflows (indicated by ellipsoidal areas). Only examples with a sub- Alfvenic disk wind 
proflle were chosen for this flgure. In general, a "flat" Alfven surface (i.e. curving towards 
the r-direction) indicates a less coUimated jet. Some of these outflows are not fully evolved 
dynamically on a global scale, but as we measure the collimation degree only within r = 150, 
the main results derived are nevertheless rehable. For highly coUimated outflows the Alfven 
surface curves parallel to the jet axis already at rather low altitudes above the disk. 
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Fig. 8. — Showcase run of Alfven surfaces for different simulation runs (see Figs. 2, 3, 6, 9). 
Selected are simulations: p20, il6; a2, ill, p4; pl9, pl6, pl5, il5; 13, 110; 12, 19; p7, 11, pl4 
(clockwise from r-axls) . The shape of the Alfven surfaces s governed by the parameters /i and 
Hp. Dashed/dotted hues highlight Alfven surfaces leaving the grid through the 2;-boundary 
/ r-boundary. As indicated by the ellipsoids, the latter ones are typical for highly coUimated 
outflows, the first ones for less coUimated outflows. Only examples with a sub-Alfvenic disk 
wind profile were chosen for this figure. 
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We have also extreme examples (i5, 18) for which the Alfven surface curves towards the 
equatorial plane (see Fig. 11 for 18) such that the outer part of the disk wind is launched 
super- Alfvenic. Those winds cannot not magneto- centrifug ally accelerated but are magneti- 
cally driven by the Lorentz force in vertical direction, Fl^z ~ jV x -B^ + x B^.. Note that 
simulations il5 and p20 seemingly look similar to the case of a super- Alfvenic disk wind, but 
actually start as sub- Alfvenic disk wind everywhere as resolved by our scaled grid. 

The kink in the Alfven surface in simulation ilO (Fig. 2, bottom) demonstrates that the 
outer part of the outflow has not yet reached a stationary state. The foot point radius of 
that flux surface intersecting the kink is at about 50 ri. Similarly, this holds for most of the 
other simulation runs. 



5. Summciry 

We have performed numerical MHD simulations of jet formation from accretion disks. 
The disk is taken as a time-independent boundary condition determining the disk magnetic 
flux proflle and the density proflle of the mass flux from the disk surface. This implies that 
we implicitly assume that disk internal physical processes such as diffusion, advection, or 
turbulence have worked together in generating and maintaining a large scale disk magnetic 
field structure which leading to a disk surface field profile and the mass flux profile as 
prescribed in our simulations. The underlying assumption is that a variation of internal disk 
parameters can lead to a variation in the resulting magnetic and mass flux profile. 

We applied the ZEUS-3D MHD code in the axisymmetric option and modified for phys- 
ical magnetic diffusivity. Our physical grid size is (150 x 300) inner disk radii corresponding 
to about (6.7 x 13.3) AU for ri ~ 10 in the case of a protostellar jet which is several 
times larger than in previous studies in the literature and comparable to the observational 
resolution. We applied a simple parameterization for physical magnetic diffusivity in which 
turbulent Alfvenic pressure is responsible for both turbulent magnetic diffusivity and pres- 
sure. However, the magnitude of magnetic diffusivity is well below the critical value derived 
earlier beyond which the degree of jet coUimation becomes affected (Fendt & Cemeljic 2002). 

The major goal of this paper was to investigate if and how the jet collimation depends on 
the accretion disk magnetic flux proflle. As a quantitative measure of the collimation degree 
we apply the directed mass flux (axial direction versus lateral direction). We have run a 
substantial number of simulations covering a large area in the parameter space concerning 
the mass flux and magnetic flux profile across the jet. In particular, we varied 

• the disk poloidal magnetic field profile parameterized by a power law Bp{r) r~^, 
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• the disk wind density profile (i.e. the density "injected" from the disk surface into the 
jet flow) parameterized by a power law Pp{r) r'^p, 

• the disk magnetic field strength and/or the total mass fiux. 

This corresponds to a variation of the mass flux profile across the jet, the magnetization 
profile across the jet, and the magnetosonic Mach number across the jet. 

Our simulations give clear evidence that a flatter magnetic fleld proflle in the disk wind 
generally leads to a higher degree of jet coUimation. This holds for all density profiles applied, 
however, the coUimation degree increases with increasing slope of the density profile. For 
very fiat magnetic field profile (// ^0.4) the jet fiows do not settle into a stationary state. 
Instead, instabilities evolve along the outfiow axis. 

Furthermore, we find a unique correlation between the slope of the disk wind magneti- 
zation profile and the degree of jet coUimation, potentially explaining why not all accretion 
disks do launch a jet and why a disk does not launch jets for all its life time. 

A variation of the total magnetic field strength does not seem to affect the degree of 
coUimation (at least not for a variation by factors up to 8 we investigated). The same holds 
for the total mass fiux in the jet fiow. However, the relative field strength affects the time 
scale for the dynamical evolution of the outfiow - jets with stronger magnetic fiux evolve 
faster. 

Comparing simulations with different magnetization but the same magnetic fiux profile 
we were able to prove the Michel scaling (Michel 1969) for the asymptotic outfiow velocity 
V(X) ~ cr^^^ with relatively good agreement. The small deviations may be explained by the 
fact that we investigate coUimated and not spherical outfiows and also outfiows which are 
to some part still dynamically evolving on the global scale. 

In order to prove the feasibihty of " poloidal coUimation" which has been discussed in 
the literature, wc invented a toy model introducing artificial decay of the toroidal magnetic 
field component (in addition to our parameterization of physical magnetic diffusivity). In 
general, such outflows show less coUimation and also reach only lower asymptotic speed. 
However, as a matter of fact, outflow coUimation is deflnitely present even in the absence of 
toroidal flelds. 

Comparing our numerical results to observations of DC Tau wc flnd indication for a 
flat magnetic fleld proflle in this source, as in this case the observed and the numerical jet 
velocity (low velocity component) and jet rotation do match. 

In summary, our paper clearly states a correlation between the disk wind magnetization 
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profile and the outflow collimation. Good coUimation requires a flat magnetization proflle 
across the jet, i.e. sufficient magnetization also for larger disk radii. The final solution to 
the question of the jet mass loading and the disk-jet magnetic field structure is expected to 
come from simulations taking into account self-consistently both the internal disk dynamics 
and the outflow. 

1 thank the LCA team and M. Norman for sharing the ZEUS code. H. Zinnecker and 
J. Wambsganss are thanked for financial support during the time in Potsdam when most of 
the computations for the present paper were performed. 
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Fig. 9. — Simulation runs with different disk boundary conditions (see Tab. 1). Density and 
poloidal magnetic field lines for il6 {t = 2040), pl5 {t = 990) il {t = 2000) p7 {t = 1000), 
al {t — 2830), from top to bottom. See Fig. 1 for further notation. 
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Fig. 10. — Simulation runs with different disk boundary conditions (see Tab. 2). Enfianced 

magnetic diffusivity for tfie toroidal magnetic field component B^p. Density and poloidal 
magnetic field lines for plO (t = 1800), pll {t = 3100), pl7 (t = 1890), il2 {t = 1520), 
pl2 {t = 870) from top to bottom. See Fig. 1 for further notation. 
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Fig. 11. — Simulation run with different disk boundary conditions (sec Tab. 1). 
{t = 290) for a parameter run witli tlie Alfven surface entering tlie disk surface, 
poloidal magnetic field lines. See Fig. 1 for further notation. 
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